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SUMMARY 

A distributed parameter model for computing unsteady flow conditions in liquid sys- 
tems is presented. The model is adapted for use in conjunction with a digital computer. 
The analytical method employed is essentially a solution synthesized from the effects of 
incremental step pressure pulses. The pressure pulses are generated because of incre- 
mental flow-rate changes that originate in a fluid system from a variety of sources, in- 
cluding the mechanical motion of the system structure. The pressure pulses propagate 
throughout the system at sonic velocity and are partially transmitted and reflected at 
each discontinuity. The velocity change caused by each pressure pulse is obtained from 
the relation for the characteristic acoustic impedance. Pressure and velocity time his- 
tories at any point in the system are obtained by a timewise summation of the contribu- 
tions of the incremental pressure pulses passing that point. 

The analysis is presented in a form general enough to be applied to a variety of 
liquid-filled fluid systems. To illustrate the application of the method to a specific sys- 
tem, the response of a long, straight propellant line to a sinusoidal inlet flow and pres- 
sure perturbation is computed. Both a constant-cross-section line and a tapered line 
are analyzed in the example. 


INTRODUCTION 

Fluid systems used in the aerospace industry often require high volume flow rates 
together with a small permissible range of variation from the design flow and pressure 
operating points. Further, the trend is towards allowing increasingly short periods for 
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the startup and shutdown transients. As fluid systems become large in size, the dy- 
namic contributions of components ordinarily not considered in the dynamic analysis of 
small systems become important. Representing liquid-filled lines, for example, by a 
single lump in a lumped-parameter analog-computer systems analysis is not adequate for 
fast transients or high-frequency sinusoidal disturbances. The calculation of unsteady 
flow containing high-frequency disturbances requires a distributed parameter model. 

Equations representing pressure and flow disturbances in a liquid network can be 
derived from continuity and momentum principles. Closed-form distributed parameter 
solutions for sinusoidal flow perturbations in long lines having negligible fluid damping 
are available and are in good agreement with experimental data (refs. 1 to 3). In gen- 
eral, however, the nonlinear partial differential equations obtained are difficult to solve 
without resorting to tedious numerical or graphical techniques. 

An analytical study was therefore undertaken at the NASA Lewis Research Center to 
develop a numerical distributed-parameter solution for unsteady flow in liquid systems 
that would be in a form readily lending itself to digital -computer computation. The ap- 
proach used to obtain the numerical solution of the continuity and momentum equations 
is fundamentally similar to that used in the method of characteristics solution (ref. 4) 
but is simpler in form. The method will be referred to herein as the wave-plan solution. 

This report presents the details of the wave-plan analysis along with the basic ele- 
ments of the corresponding digital-computer program required to calculate unsteady flow 
in liquid systems. The analysis is presented in a form general enough to be applied to a 
variety of fluid systems. Examples are given to illustrate the application of the method 
to specific liquid systems. 


ANALYSIS 

General 

Basic equations . - The pressure and velocity of the liquid within a line as a function 
of position and time can be obtained from the equations of momentum and continuity. 

The momentum equation for a one -dimensional elastic fluid of constant mean density 
is (ref. 4) 
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( 1 ) 


where /(v) is the resistance due to fluid viscosity, which is some function of velocity. 
(Symbols are defined in appendix D. ) 
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The continuity equation for an elastic liquid in an elastic line is (ref. 4) 


3H + y 3H _ .C 2 3V 
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Equations (1) and (2) are nonlinear partial differential equations which, when solved 
along with the boundary conditions in a liquid-filled line, give relations for velocity and 
pressure in the line. 

Method of solution . - The direct solution of the continuity and momentum equations 
to obtain velocities and pressures is difficult even when the boundary conditions are rela- 
tively simple (open or closed ends). In addition, exact solutions have been found only for 
a system where the resistance term is linearized or neglected entirely (e. g. , refs. 1 
to 3). 

The continuity and momentum equations can be solved implicitly by employing a 
numerical process similar to a method of characteristics solution. This method, re- 
ferred to as a wave -plan solution, provides a straightforward, accurate distributed 
parameter solution for unsteady flow in a liquid system. 

A wave -plan solution is obtained as follows: At the point in a liquid system where a 
disturbance is introduced (such as an oscillating valve or a moving impedance change), 
an incremental change in liquid flow rate due to the disturbance over a very short inter- 
val of time is computed. The incremental pressure pulse accompanying the flow-rate 
change is then computed. This pressure pulse is propagated throughout the system at 
sonic speed. The pressure pulse is partially reflected and partially transmitted at all 
geometrical and physical discontinuities in the fluid network. Pressure and velocity 
time histories are computed for any point in the system by summing with time the con- 
tributions of incremental waves. 

Characteristic impedance . - The relation (characteristic impedance) between pres- 
sure and velocity changes caused by a pulse traveling in the liquid-filled line is computed 
from momentum considerations. Figure 1 shows pressure and flow conditions, as a 
pressure wave propagates in a liquid-filled line, that exist at time t and at time t + At. 

The wave takes the time At to travel the dis- 
tance Ax from A to B. During this time there 
is a pressure P + AP on the left end and a 
pressure P on the right end of the liquid con- 
tained between points A and B. This unbalanced 
pressure causes the fluid to accelerate. 

Newton's second law gives 

(P + AP - P)A = pAAx^ 

At 
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Figure 1. - Effect of pressure pulse on mean line conditions. 
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Canceling and rearranging yield 


AP = p AV — 

At 

The term Ax/ At is the propagation speed of the pressure wave. The wave speed is 
equal to the sonic velocity C in the system if the mean velocity of the liquid in the line 
is neglected. Since the mean velocity of the liquid is usually much smaller than the 
acoustic velocity, this is usually permissible. Thus, 

AP = pC AV (3a) 


or in terms of head of liquid 


AH = — AV (3b) 

g 

The sonic speed C for a liquid flowing within a line is influenced by the elasticity of the 
line wall, and for a system that is axially unrestrained (refs. 4 to 6) can be calculated 
from 


C = 



Equation (3), which is the well known Joukowski equation, states that the pressure 
change at any point in a line is equal to the product of the velocity change at the point 
times the characteristic acoustic impedance pC of the liquid in the line. This relation 
may also be derived by employing continuity and energy principles (refs. 5 and 6). 


Generation of Pressure Waves 

Pressure waves may originate in liquid flow systems in a variety of ways. Two 
common sources are an orifice or valve with a varying open area and a moving point in a 
line where there is an impedance change. In order to apply the wave-plan method to 
varying area change or moving impedance change, the corresponding disturbance func- 
tion is approximated by a series of discrete changes over small equal time intervals. 
Figure 2 shows this approximation. 

Each of these changes is assumed to occur instantaneously at some time during the 
time interval. The pressure perturbations associated with each of these changes 
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are then computed as shown in the following examples. 

Varying area orifice . - A pressure wave is gen- 
erated at an orifice that experiences a change in open- 
ing area. The magnitude of the generated wave is a 
function of conditions before the area change and the 
magnitude of the area change. 

The flow through an orifice is assumed to obey the 
square -law relation: 

v = b Yah 

The orifice coefficient B is a function of the open area, the line area, and the dis- 
charge coefficient and is easily determined for a particular orifice. For the case of an 
open area that is a small percentage of the line area (Aq « A), the orifice coefficient is 
given closely by 

where the discharge coefficient is a function of Reynolds Number (and hence, veloc- 
ity). For small velocity perturbations, however, the discharge coefficient varies only 
slightly, and hence the orifice coefficient may be considered to vary linearly with area 
changes. The examples in this report deal with this case. 

Figure 3 shows conditions at an orifice before and after a small orifice-area change. 
The flow is from left to right. The external pressure head may be varied in a pre- 
scribed manner. 

The momentum equation across the wave front gives 

iH U = ^< V l - v u> 



Time, t 


Figure 2. - Step approximation of disturbing 
function. 


The flow out the orifice after the orifice area change is given by 


1 Time 

Vj, Hj |h 2 t 

B 1 

H AH U 

Vj, Hj V n , H u ^ t + At 

J B 2 

Figure 3. - Effect of step change in orifice coeffi- 
cient at terminal orifice. 


V 11 - B 2 V H 1 + AH U - H 22 

Solving the momentum and orifice equations simulta- 
neously gives the following quadratic in Vj^ 

V 2 n + bV n + c = 0 (5) 

where 
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o / CV 
c = -b|(h 1 + — 

The positive root of this equation gives the desired value for outflow velocity. Substi- 
tuting this value for Vjj back into the momentum equation gives the magnitude of the 
pressure wave generated by a sudden change in orifice coefficient from to Bg. 

Moving terminal orifice . - The analysis of the magnitude of a pressure wave gen- 
erated by the motion of a component in a fluid line is exemplified by considering the 
motion of a terminal orifice. The lateral velocity of the orifice is so represented by a 
series of step changes that the velocity of the orifice remains constant over short time 
intervals. Figure 4 shows conditions at a terminal orifice just before and a short time 
after a step change in velocity. The orifice coefficient need not be constant. 

The momentum equation across the generated wave is 

iHn=|(Vi-V u ) 



The displacement shown in figure 4 takes place in the time increment At. The con- 
tinuity equation states that the net inflow across boundary AA equals the net flow out the 
orifice plus the storage that takes place during the time interval. Therefore 

V U A At = VE 2 A At + VO A At 


or 


V 11 = VE 2 + VO 
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where VE 2 is the lateral velocity of the ori- 
fice end during the time interval, and VO is 
the liquid velocity in the line (with respect to 
the orifice). 

The flow out the orifice is given by 
VO = B 2 ^/Hj + AH n - H 22 


Figure 4. - Effect of lateral motion of terminal orifice. Solving the preceding equations simultaneously 
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for the resulting line velocity gives the following quadratic in VO: 


VO 2 + bVO + c = 0 


( 6 ) 


where 

BpC 
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Effect of Viscous Resistance 

After the generation of a pressure wave at a perturbing point in a liquid system, the 
wave propagates with sonic velocity throughout the system. The viscous resistance of 
the liquid medium influences the propagation of this wave. 

The effect of viscous resistance on pressure pulse propagation can be neglected 
without appreciable effect for short, large-diameter lines where such losses are small; 
however, for longer lines or small-diameter lines carrying viscous liquid, the resistive 
losses are not negligible and should be included in an unsteady flow analysis. 

The exact solution of the partial differential equations (1) and (2) is limited by the 
necessity of linearizing or neglecting the friction term entirely. For graphical analyses 
it has been necessary to consider the friction losses as lumped at one or more points in 
the pipeline if a manageable solution is to be obtained (ref. 7). 

To include the effects of viscous resistance in a wave-plan solution, the extent that 
fluid viscosity influences the propagation of a pressure pulse must be determined. 

There is presently no experimental or analytical information available which is in a 
form that can be used for the prediction of the effect of viscous resistance in unsteady 
line flow. Because of this limitation it is necessary to make a quasi-steady approxima- 
tion and to assume that viscous losses in a small line length dx are given at any instant 
by the Darcy equation as 


Ah L = 


f Ax V 2 
2gD 


(7) 


The form of equation (7) is identical to that of an internal square-law orifice if the 
orifice coefficient is given by 
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This equation implies that viscous losses over a small length of line can be lumped 
at an internal square-law orifice with a properly chosen orifice coefficient. This repre- 
sentation is referred to as the “orifice analogy, ** and its use in graphical analysis has 
been suggested by Bergeron (ref. 7). 

Line losses can be distributed at many discrete points along a line by inserting a 
large number of correctly chosen square -law (friction) orifices. Impinging waves are 
then reflected and transmitted at these friction orifices in a manner similar to that of 
reflection and transmission through a small region of flowing viscous liquid. The equa- 
tions governing wave reflection and transmission at a friction orifice are developed in 
the next section. 


Reflection and Transmission of Pressure Waves at System Discontinuities 


It is necessary at each discontinuity to compute the magnitude of transmission and 
reflection of each pressure pulse in terms of initial conditions at a discontinuity, the 
nature of the discontinuity, and the magnitude of the impinging waves. These deriva- 
tions are all made by employing the three fundamental fluid flow relations (continuity, 
momentum, and energy). 

Figure 5 shows two waves impinging and reflecting simultaneously at a discontinuity 
in a liquid system. The notation used in this figure is employed in subsequent specific 
examples. Liquid moving from left to right has a positive magnitude of velocity. 

The momentum equations across the wave fronts yield 



Time 
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t + 


Figure 5. - Nomenclature for conditions before and after wave 
action at discontinuity. 


(See fig. 5 for definition of V* and V M .) 
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AH 22 = X (V 22 - V "> 

Equations (8) and (9) combine to give 

+ Y (v i " v n ) 
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and equations (10) and (11) give 

4H 22 = ah 2 + -2 (v 22 - v 2 ) 
s 

The new pressures are given by 


H 11 = H l + AH 1 + AH 11 


(13) 


H 22 “ H 2 + AH 2 + AH 22 


(14) 

(15) 


Terminal orifice . - A terminal orifice is considered to be bounded by a pressure 
reservoir. The pressure in the reservoir may be changing in a prescribed manner. In 
addition, the orifice coefficient may also be changing in a prescribed manner. Figure 6 
shows conditions before and after reflection from an orifice bounded at the inlet by a 
pressure reservoir. 

If the discharge after wave action is in the posi- 
tive direction (from the reservoir to the line), the 
head-discharge relation for the orifice is given by 
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Figure 6 . - Pressure-pulse reflection at terminal 
orifice. 


V 22 = b 2 iSl - < h 2 + ah 2 + aH 22> 

Solving this equation simultaneously with the 
momentum equations (13) for the velocity gives the 
following quadratic in V 22' 


V 22 + bV 22 + c = 0 


(16) 


where 


b = 


b|c 2 

g 


and 


c = - 


2 / C qV o' 

B 2 (hi “ H 2 “ 2 ah 2 + 


Because the resulting direction of flow was taken as positive, the positive root of 
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equation (16) is the desired result. 

The magnitude of the reflected pressure wave is obtained from equation (13). The 
new pressure at the orifice is given by equation (15). 

If the direction of flow is into the reservoir after wave action (negative direction), 
the velocity head relation is given by 

V 22 = B 2 V<H 2 + AH 2 + AH 22 ) - H n 


Simultaneous solution with the momentum equations gives a quadratic in V 22 : 

V22 + W 22 +c = 0 (17) 

where 

h _ B I C 2 

g 


c 




- 2 AH 2 



The negative root of equation (17) is the desired result. 

It is necessary to determine the resulting direction of flow so that the correct equa- 
tion (eq. (16) or (17)) can be applied. Inspection of equation (16) discloses that the term 
Hu - H 2 - 2 AH, + C 2 V 2 /g must be positive to yield the necessary positive root. If 
this term is negative, the flow must be into the reservoir, equation (17) must be used, 
and the negative root must be computed. 

Friction orifice . - Figure 7 shows conditions before and after pulse-wave action at 
a friction orifice. Because of identical conduits on both sides of the orifice, the head- 
velocity relation for the orifice after wave action for flow from left to right is 

V 11 - v 22 = Bf V H 1 + AH 1 + iH U - H 2 - aH 2 - iH 22 



r^AH! — 
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Figure 7. - Pressure-pulse reflection at internal orifice. 


Noting that = C 2 = C and solving this equation 
simultaneously with the momentum equations (12) 
and (13) give 


V 11 + bV ll + c = 0 


(18) 


where 
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b = 


2B|c 
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c = 



Hj + 2 AHj - 


Hr 


- 2 AH 2 + 



The positive root of this quadratic equation gives the resulting velocity through the 
friction orifice. The term + 2 AH^ - H 2 - 2 AH 2 + 2CVj/g must be positive to yield 
a positive root. If this term is negative, the resulting flow is from right to left, and the 
head-velocity relation for the orifice is 

V 11 = V 22 = B F V H 2 + AH 2 + AH 22 " H 1 " AH 1 “ AH 11 


This equation combines with the momentum equations to give the following: 

V 11 + bv H + c = 0 


(19) 


where the coefficients b and c are of the same magnitude and of opposite sign from 
the corresponding coefficients in equation (18). 

The negative root of this equation gives the resulting velocity in the negative direc- 
tion. The magnitudes of the resulting pressure waves and the pressures after wave 
action are given by equations (12) to (15). 

Diameter discontinuity . - The relation between net head and discharge for a diam- 
eter change in a conduit is derived by utilizing energy relations. Figure 8 shows the 
notation for this derivation. 

If Aj > A 2 , the diameter discontinuity is an abrupt constriction, and the energy 
equation for flow through the constriction from left to right after wave action is 
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Figure 8. - Pressure-pulse reflection at diameter discontinuity. 
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The area ratio is denoted as 

A 1 
R = — - 


Then the continuity equation is 


^2\^k 

aJ 2g 


( 20 ) 
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( 21 ) 


V 22 = R V n 


If these relations are substituted into the energy equation and the terms are com- 
bined, the following net head-velocity relation is obtained: 


o / 3R - 2R - 1 
H H ~ H 22 = ^lly~ 4g j 


( 22 ) 


If A^ < Ag, the diameter discontinuity is an abrupt expansion. The energy equation 
for flow through the expansion from left to right after wave action is 


111 

2g 
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This equation gives the following velocity-head relation for an expansion: 
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(23) 


The case of a lossless diameter discontinuity is also important because step diam- 
eter changes can be employed to simulate a tapered line. In the continuously tapered 
line the losses are small and sometimes can be neglected. The energy equation without 
the loss term is 



+ H 11 = 
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22 


2g 
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22 


The resulting velocity-head relation is 


H 11 ' H 22 = V 11 


2 /R" 


2g 


The final net head after wave action can be written as 


H H " H 22 = % + AH X + AH n - H, - AH, - AH 22 
or 


( 24 ) 
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( 25 ) 




H 11 “ H 22 = H 1 + 2 AH 1 + — (V t - V n ) - H 2 - 2 AH 2 - ^ (V22 _ yj 

g g 

The relation is combined with the continuity equation and the velocity-head relation 
to obtain second-order polynomials for the line velocity after wave action. The equation 
is 


aV ll + bV ll + c = 0 


(26) 


where 


C. C 0 R 
g g 
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The coefficient a depends on the type of discontinuity and whether energy losses 
are included. For a contraction with losses included the coefficient is 


3IT - 2R - 1 
4g 


(27) 


For an expansion with losses included the coefficient is 


a = 


R 2 - R 
g 


(28) 


For either type of discontinuity with no losses included the coefficient is 


a = 


R 2 - 1 

2g 


(29) 


For each case the resulting velocity is the positive root of equation (26). There- 
fore, the term 


C 1 V 1 C 2 V 2 

Hj + 2 AH. + — — - + — ^ - H 2 - 2 AH 2 
g g 
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must be positive. If this term is negative, the resulting flow is from right to left, and 
the equations are reformulated by putting the loss term on the correct side of the energy 
equation. 

Junctions of thr ee or more legs. - The previously presented mathematical analyses 
pertaining to the reflection of a pressure pulse at a diameter discontinuity have taken 
into account energy losses introduced by the discontinuity. With junctions of three or 
more legs, the inclusion of energy loss terms becomes increasingly difficult because of 
the irreversibility of frictional effects. Thus, the treatment employed herein will be 
similar to the one presented in standard references on fluid dynamics and will not be 
derived (e. g. , ref. 8). 

The relations for computing the percentages of magnitude of a pressure wave trans- 
mitted and reflected at a junction are based on the following: 

(1) The Joukowski equation applied across each wave 

(2) Continuity of flow at the junction 

(3) Continuity of pressure at the junction 

The magnitude of the wave that is transmitted to all the other legs of a junction of 
n legs due to a wave of magnitude AH impinging in line i is given by T(i)AH, where 
the transmission coefficient T(i) is given by 


T(i) = 


2A(i) 

C(i) 



m 

C(j) 


(30) 


The magnitude of the wave reflected in leg i is given by R(i)AH when the reflec- 
tion coefficient R(i) is given by 


R(i) = T(i) - 1 

Analysis of the Dynamics of a Liquid System by Wave Plan 

The wave-plan analysis supposes a system composed of a discrete number of dis- 
continuities connected by lossless line segments, which serve only to transmit pressure 
pulses. The discontinuities include geometric ones, such as terminal orifices and diam- 
eter discontinuities, and artificial ones, such as friction orifices. A typical simulated 
fluid system is shown in figure 9. 
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Figure 9. - Wave-plan representation of typical fluid system. 

There is a variable-area orifice at A (valve) and a fixed orifice at B. At C there is 
an abrupt change in line diameter. Friction orifices are inserted between A and C and 
between C and B to simulate viscous resistance. 

The equations developed thus far have applied only to conditions at a particular point 
in a system at a particular instant of time. Through the use of these equations it is 
possible to compute velocity, pressures, and magnitudes of pressure pulses leaving a 
particular discontinuity in terms of the magnitudes of impinging waves and conditions 
prior to wave action. 

To analyze a system, however, the equations must be both time and position de- 
pendent. This dependency is indicated by the form of the basic partial differential equa- 
tions (eqs. (1) and (2)). It is apparent that the pressure wave impinging at one discon- 
tinuity in a liquid system is exactly that wave which left an adjacent discontinuity at some 
time in the past (because of the lossless line connector simulation). Thus, the waves 
are related by position and time. In general, variables have the form /(x, t). Specifi- 
cally, H(x,t), V(x,t), and AH(x, t) denote the pressure head, velocity, and pressure 
pulse as a function of position and time, with x and t being the position and time sub- 
scripts, respectively. 

The wave -plan analysis provides for information at specified points in the system 
(discontinuities) and at discrete time intervals. New information is available only at 
discrete time intervals because all of the disturbing functions are approximated by a 
series of small step changes occurring at specified times. The system response to 
these disturbing functions will also be in the form of step changes. Because of the two- 
parameter dependence of the dependent variables, it is advantageous to denote position 
and time subscripts. 

Position subscripts . - In order to apply the analysis equations that have been devel- 
oped for a discrete point to a liquid system, it is necessary to introduce a method of 
denoting position. This is done by numbering adjacent discontinuities in a consecutive 
manner as indicated in figure 9. (This arrangement is modified at junctions of three or 
more legs). The position subscript is denoted as L. 

In addition, because pressure and velocity states may differ across discontinuities, 
it is necessary to denote a left and a right side of a discontinuity. This is done by adding 
an L for left and an R for right to the nonsubscripted part of the variable. 

Thus HL(L, t) denotes the pressure head to the left of discontinuity L at time t. 
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Time subscripts. - The time subscript J is defined as follows: 


t = J At 

where t is the time and At is the working time increment. Thus HR(L, J) denotes the 
pressure head to the right of discontinuity L at time t = J At. 

Selection of w orking time increment. - The working time increment At is the time 
interval between successive computations. 

Its selection is determined by two factors. First, the time increment must be 
small enough to approximate accurately all disturbing functions by a series of step 
changes. Second, it is necessary that all reflections of pressure pulses in the system 
take place at ah integer number of time intervals. The wave travel times between all 
adjacent discontinuities will then be an integer number of working time increments. If 
this were not the case, waves would be impinging on a discontinuity in a completely 
arbitrary fashion and would make the solution unmanageable. 

The selection of the time increment is simplified because the majority of discon- 
tinuities in a system will be friction orifices, and these can be placed where desired. 

For example, the selection of the working time interval for the system shown in figure 9 
is made as follows: 

(1) The wave travel times between adjacent geometric discontinuities are determined: 

(a) tac - wave travel time between A and C 

(b) tcb - wave travel time between C and B 

Within the desired limits, the largest time interval of which these travel times are 
integer multiples is determined. This time interval represents the largest working time 
increment possible. 

(2) It is then determined if the integer number of working time increments necessary 
for travel between A and C and between C and B is large enough to insert the desired 
number of friction orifices that must, of course, be separated in time by at least one 
working time increment. If not, the increment can be divided by any integer to obtain a 
smaller working time increment. 

(3) It must be determined if the working time interval necessary to assure an integer 
number of time intervals between all discontinuities is small enough to approximate 
accurately the disturbing function (variable -area orifice at A) by a series of step changes. 
If not, the working time increment is further divided by an integer necessary to get a 
suitable approximation. The integer number of time increments between discontinuities 
is increased by this factor. 

Subscripted n ota tion for analysis equ ations. - The analysis equations are easily ex- 
tended for application to a liquid system by the use of the subscripted notation. Fig- 
ure 10 shows conditions at a discontinuity before and after wave action when the sub- 
scripted notation is employed. 
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Figure 10. - Subscripted nomenclature for conditions before and after wave action at 
discontinuity. 


The time that the waves reach discontinuity L is t = J At. Conditions before and 
after that time are constant, and step changes occur at t = J At. 

A comparison of the nonsubscripted notation in figure 5 to the subscripted notation 
in figure 10 gives the identities necessary to make the analysis equations general expres- 
sions for a liquid system. These identities are as follows: 

VL(L, J - 1) = = velocity to left of discontinuity L 

HL(L, J - 1) = H-^ = pressure head to left of discontinuity L 
VR(L, J - 1) = V 2 = velocity to right of discontinuity L 
HR(L, J - 1) = H£ = pressure head to right of discontinuity L 

AHR(L - 1, J - KX) = AHj = pressure pulse coming from adjacent discontinuity to 
left. KX is the number of working time intervals it takes a sonic disturbance 
to travel between the two discontinuities. 

AHL(L + 1, J - KY) = AH 2 = pressure pulse coming from adjacent discontinuity to 
right. KY is the number of working time intervals it takes a sonic disturbance 
to travel between these two discontinuities. 

The conditions after wave action are as follows: 

VL(L, J) = Vjj = velocity to left of discontinuity L 
HL(L, J) = H n = pressure head to left of discontinuity L 
VR(L, J) = V 22 = velocity to right of discontinuity L 
HR(L, J) = H 22 = pressure head to right of discontinuity L 

AHL(L, J) = AH-q = pressure pulse leaving discontinuity L and moving to left 
AHR(L, J) = AH 22 = pressure pulse leaving discontinuity L and moving to right 

Substitution of these identities into the analysis equations (eqs. (8) to (29)) yields a 
perfectly general set of algebraic equations applicable to any liquid system. 

Method of computation . - The solution is carried out by using the subscripted equa- 
tions as follows: 

(1) A working time increment is selected. 

(2) Initial conditions are computed from given steady- state values that exist in the 
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system prior to the initiation of a disturbance. These conditions include velocity and 
pressure to the left and right of each discontinuity and are denoted at t = 0 by VL(L,0), 
VR(L,0), HL(L,0), and HR(L,0). 

(3) Computations are then carried out to determine conditions at each discontinuity 
at the end of the first working time interval by using the conditions from step (2) as ini- 
tial conditions. It should be noted that the analysis equations will predict no change at an 
undisturbed discontinuity that is not subjected to impinging waves. Thus, conditions will 
only change at the end of the first time interval at a disturbance source (point A, fig. 9). 
At the disturbance source conditions change, and waves are generated that will subse- 
quently alter conditions at adjacent discontinuities. 

(4) Each discontinuity is analyzed at the end of each interval. The computations are 
carried out as long as desired to give pressures and velocities at each discontinuity at 
the end of each working time interval. 


Digital Computer Programing 

The wave-plan analysis consists of the sequential solution of many equations. The 
dynamic analysis of even a rather simple liquid system for a reasonable number of time 
intervals involves a large number of computations. A digital computer must be employed 
if the solution is to be obtained in a reasonable length of time. 

To make the digital analysis as general as possible, computer segment programs 
(subroutines) have been written for the most common types of discontinuities. These 
discontinuity solutions can then be combined (incorporating junction equations, if neces- 
sary) to obtain digital computer models for various liquid flow systems. 

Each discontinuity subroutine computes conditions at the discontinuity for some 
point in time as a function of the local velocity and pressure-head conditions a time in- 
terval earlier, the magnitude of impinging waves (from adjacent discontinuities), and the 
step changes in disturbing functions during the time interval. 

The following subroutines have been written to be used in the digital model repre- 
sentation: 

(1) Terminal orifice subroutine 

(2) Internal (friction) orifice subroutine 

(3) Diameter discontinuity subroutine 

Each subroutine uses the general relations given by equations (12) to (15). Also 
each subroutine initially assigns nonsubscripted identities to their subscripted counter- 
parts and finally assigns the nonsubscripted results to the subscripted terms (as indi- 
cated in the section entitled Subscripted notation for analysis equations). A description 
and computer printouts of the digital computer subroutines are included in appendix A. 
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The subroutines are combined to form digital computer models for the dynamic 
analysis of liquid flow systems. Two examples are included in the next section. 

RESULTS AND DISCUSSION 

The equations and programs developed in the ANALYSIS section form a basis for the 
analysis of unsteady flow conditions in a liquid flow network. The equations and pro- 
grams have been kept general so that they may be applied to the analysis of a variety of 
liquid systems. In general, the wave -plan digital- computer analysis has distinct advan- 
tages over other methods such as analog- computer models that substantially linearize or 
lump parameters and small perturbation techniques that linearize around mean-line con- 
ditions. The following features of the wave-plan analysis should be noted: 

(1) Viscous friction effects are easily included. 

(2) Perturbing functions of any form may be used. 

(3) Nonlinear relations are easily included. (Linearization of parameters is of little 
or no advantage. ) 

(4) Both transient and steady-state responses to disturbing functions are given. 

(5) Complex networks can be analyzed. 

Two examples will be used to illustrate the application of the wave-plan method to 
the construction of digital models of particular liquid flow systems. The first example 
will illustrate the method of computing the transient and steady-state response of a long, 
straight hydraulic line to a periodic sinusoidal flow disturbance. In the second example 
the response of a long tapered line will be computed. 


Example 1: Long Hydraulic Line With Oscillating Inlet Orifice 

This example was chosen because the computer results can be compared directly 
with experimental data from the Lewis line dynamics rig (refs. 1 to 3). Schematics of 
the liquid flow system are given in figure 11. Figures 11(a) and (b) show the analytical 
model and the experimental rig, respectively. The system consists of a 68-foot-long, 
7/8-inch-inside-diameter line with pressure reservoirs at both ends. The test fluid was 
JP-4 fuel or another hydrocarbon, depending on the liquid viscosity desired. Line pres- 
sures ranged from 200 to 400 pounds per square inch. A variable-area throttle valve 
that perturbs the system is located at the upstream end of the line (point A), and a fixed 
orifice is located at the downstream end (point B). The details of constructing the digi- 
tal model are given in appendix B along with the program constants used in the calcula- 
tions. 
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Figure 11. - Hydraulic line with oscillating inlet orifice. 


Effect of amplitude o f input perturba tions. - Most line-dynamics analyses are based 
on small-perturbation techniques that permit the use of linearized terminal impedances 
as an approximation for nonlinear pressure-flow relations at the end of the line. There- 
fore, in taking experimental line -dynamics data (refs. 1 to 3) the practice is to limit the 
amplitude of the oscillating throttle disturbance generator. The minimum usable ampli- 
tude, however, for a given line condition and terminal impedance is limited in practice 
by the ratio of the sine wave signal to the nonharmonic noise which decreases with de- 
creasing amplitude. The permissible maximum amplitude is limited in practice by the 
appearance of harmonics in the sinusoidal pressure signal as the amplitude is increased. 
The system must be operated between these two limits. 

The relation between the disturbance amplitude and the degree of nonlinearity of the 
pressure and velocity perturbations was determined analytically by varying the amplitude 
of the sinusoidal input perturbations for the analytical line model. The mean orifice 
coefficients were 0. 8 at the inlet and 1. 0 at the outlet. The steady line velocity was 
14 feet per second. For the digital computer program (appendix B) the amplitude of the 
input orifice coefficient perturbation BA was varied at 40 cps as shown in table I. 
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The calculated steady-state responses for the 
pressure and fluid velocity at the inlet and outlet of 
the line (points A and B, fig. 11(a)) are shown in 
figure 12 (p. 22). The velocities and pressures are 
dimensionless, having been divided by steady line 
values. Figure 12 shows that, as the amplitude of 
the disturbing function is decreased, the system 
responses are more nearly described by a 
sinusoidal -type periodic function; however, for the 
particular line condition and terminal impedance 
used in the example, nonlinearities are evident 
even for relatively small perturbations. For ex- 
ample, the inlet velocity perturbation is very clearly nonsinusoidal for an 18. 75-percent 
disturbing amplitude even though the perturbation amplitude is only 2 percent in the posi- 
tive direction and 4 percent in the negative direction. 

Transi ent response . - Dynamic systems analyses are usually based on steady-state 
responses to steady sinusoidal inputs. This analysis technique does not, however, pre- 
dict the transient response to the initiation of a disturbance or a change in the disturb- 
ance characteristics in a liquid flow system. The wave-plan solution gives both trans- 
ient and steady- state solutions. This point is clarified by examining both the transient 
and final steady-state responses to the initiation of a sinusoidal disturbing function in a 
liquid system under steady mean-flow conditions. A value for BA of 0. 15 (case 4, table I) 
was chosen for the amplitude of the disturbance function. The line response starting with 
the initiation of the disturbing function is given in figure 13 (pp. 23 and 24). For the case 
studied several cycles were required for the perturbation velocities and pressures to ap- 
proach the steady-state values. In some cases the amplitudes of the transient responses 
are considerably different from the steady responses. For example, during the early 
part of the transient period the magnitude of the inlet velocity perturbation was nearly 
double the steady-state value. 

Comparison of analytical and exp erimental results . - The experimental hydraulic 
line (fig. 11(b)) was operated with the upstream throttle valve sinusoidally perturbed at 
70 cps with area perturbation amplitudes of 12. 3, 33. 8, 55. 3, and 76. 8 percent of the 
open area. The mean value of the coefficient for the upstream orifice was 0. 65, and the 
corresponding value for the downstream orifice was 0. 42. The downstream orifice was 
slightly open with respect to characteristic. Figure 14 (p. 25) shows the resulting ex- 
perimental pressure-time oscilloscope traces for the quarter-length point of the line 
(measured from the upstream throttle valve). 

The experimental conditions were duplicated by the digital model, and figure 15 
(p. 25) shows the theoretical results for these four different input amplitudes. The 


TABLE I. - AMPLITUDES OF ORIFICE 
COEFFICIENT PERTURBATIONS 


Case 

Input orifice 

Percent of steady- 


coefficient 

state orifice 


perturbation 

amplitude, 

BA 

coefficient 

1 

0. 48 

60 

2 

.32 

40 

3 

.20 

25 

4 

.15 

18.75 
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Time increment 

Figure 12. - Effect of amplitude of inlet orifice area perturbation on steady-state response of long line. 


agreement between the wave shapes determined analytically and experimentally is very 
good for all cases. In addition the analytical and experimental magnitudes agree to 
within the accuracy of the measurements (approx. 3 percent). 

This example illustrates that the wave -plan analysis for unsteady liquid flow is 
capable of accurately including nonlinear effects and giving results in good agreement 
with experimentally observed values. 
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Time increment 

Figure 13. - Transient response of long line to initiation of sinusoidal inlet valve area perturbation. 

Example 2: Response of Long Tapered Line to Sinusoidal Input Disturbance 

The analysis of the response of a liquid system composed of lines of nonuniform 
diameter (tapered) to a periodic disturbing function is difficult if classical methods are 
employed; however, tapered lines can easily be included in a wave-plan analysis by 
approximating the tapered line by a line composed of a discrete number of diameter 
changes (fig. 16, p. 26). 

In order to evaluate the quantitative effect of tapered lines, the system shown in 
figure 17 (p. 26) was studied in detail. This system consists of two pressure reservoirs 
connected by a tapered line. At the end of the tapered lines, relatively large resistive 
losses were introduced in the form of square-law orifices. The throttle valve at point A 
has a sinusoidally varying area, and the orifice at point B is fixed. 
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Pressure 


7 Ib/sq in. 


7 Ib/sq in.- 




(b) Orifice area perturbation ampli- 
tude, 33.8 percent of open area; 
maximum amplitude (peak to peak) 
of pressure perturbation, 32.2 ± 
0.7 pounds per square inch. 


35 Ib/sq in.- 


(a) Orifice area perturbation ampli- 
tude, 12.3 percent of open area; 
maximum amplitude (peak to peak) 
of pressure perturbation, 11.2 ± 
0.7 pounds per square inch. 


35 Ib/sq in. 



(c) Orifice area perturbation ampli- 
tude, 55.3 percent of open area; 
maximum amplitude (peak to peak) 
of pressure perturbation, 73.5± 
3.5 pounds per square inch. 



(d) Orifice area perturbation am- 
plitude, 76.8 percent of open 
area; maximum amplitude (peak 
to peak) of pressure perturbation, 
147 ±3. 5 pounds per square inch. 


Figure 14. - Experimental pressure-time oscilloscope traces for quarter-length point of line. 
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(a) Orifice area perturbation amplitude, 12.3 percent of open 
area. 
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(c) Orifice area perturbation amplitude, 55.3 percent of open 
area. 
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(b) Orifice area perturbation amplitude, 33. 8 percent of open (d) Orifice area perturbation amplitude, 76. 8 percent of open 
area. area. 

Figure 15. - Analytical pressure-time traces for quarter-length point of line. 









Figure 16. - Step-diameter-change approximation to tapered line. 


This system was studied for various 
degrees of taper. The configurations 
were chosen to give the diameter at the 
inlet as 

DA = D - 6 ft (31a) 


and the diameter at the outlet as 


DB = D + 6 ft (31b) 


The diameter of the tapered conduit 
is assumed to vary linearly between the 
ends, and the average diameter is always 
D feet. 

If the resistances of the terminal orifices are large compared to the line loss, the 
steady discharge of the system is practically independent of 6 and would depend mainly 
on the reservoir pressures and the mean line diameter. Even for lines with much 
smaller end-resistive losses the system discharge would be only slightly dependent on 
the degree of taper, because line losses consisting of friction and expansion or contrac- 
tion losses would be mainly dependent on line length and average diameter. Mathemati- 
cally this can be written as 



N 

Point A B 

Figure 17. - Reservoir conduit system with tapered line. 


q 2 _ HEN - HEX 
" K en + K l + K ex 


(32) 


where K , K,, and K are the loss coefficients of the entrance orifice, the line, and 
6n a 0 x 

the exit orifice, respectively. 

When this system is used with large values of K 0n and K gx as compared to Kp 
it is possible to compare directly the dynamic response of systems that are essentially 
statically equivalent. Pressure and flow perturbations can be quantitively compared as 
a function of degree of taper. The straight conduit (6 = 0) is a special case that is in- 
cluded in the analyses. 

The tapered conduit was chosen so that the wave velocity would be constant; thus the 
ratio of diameter to wall thickness remained constant along the length of the line. 

For the digital model the tapered line was approximated by N - 1 short straight 
sections of conduit, as shown in figure 17. In this manner the diameter of each approxi- 
mating section I is computed by the equation 
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The details of the digital model along with the computer program and program con- 
stants are given in appendix C. One basic system was studied in detail, a system having 
a high-resistance orifice (three-fourths to seven-eighths closed geometrically but 
slightly open with respect to characteristic impedance) at the inlet with a 10-percent 
sinusoidal variation in the area of the orifice opening. The outlet orifice was chosen of 
even higher resistance than the inlet and represents an end that is closed with respect to 
characteristic impedance. Nine degrees of taper were studied, the average diameter D 
being 1 foot for each case. The values of 5 that were used were -0. 3, -0. 2, -0. 1, 0, 

0. 1, 0. 2, 0. 3, 0. 4, and 0. 5 foot. 

The quarter-wave resonance frequency for the nontapered line is 15 cps. Frequen- 
cies of 5, 7. 5, 10, 12. 5, 15, 17. 5, 20, and 25 cps were studied for each case. In all 
cases the pressure and flow perturbations were slightly nonlinear even though they were 
small. In all cases they were within 2 percent of mean-line values. Figure 18 shows 
typical pressure and flow perturbations for one cycle. 
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taper. 

Figure 19. - Frequency dependence of tapered-line dynamic response. 

In spite of the nonlinearity of the perturbations, they were steady and repeatable 
after several cycles (during which the transient died out) and resembled sinusoids. It 
was considered appropriate to make comparisons based on positive amplitudes. This 
procedure is followed in the ensuing discussion. 

The analytical results are displayed graphically in figure 19. Figure 19(a) shows 
the variation of the amplitude of the outlet pressure perturbation with frequency for the 
expanding taper (line that expands from inlet to outlet). Figure 19(b) shows the pressure 
perturbation gain (ratio of outlet pressure perturbation amplitude to inlet pressure per- 
turbation amplitude) for the expanding taper. Figures 19(c) and (d) show the same re- 
sults for contracting tapers. The most striking result apparent in figure 19 is the shift 
in the resonant frequency for tapered lines. This shift is shown in figure 20. 

Another result is that the maximum pressure perturbations at the outlet have a 
strong dependence on 6 and vary considerably over the entire frequency range. Thus a 
certain amount of control can be exerted over the response of this system by the judi- 
cious choice of tapers. For example, in the case of an expanding taper of 6 = 0. 3 the 
amplitude of the outlet pressure perturbation is significantly reduced when compared to 
a straight-line response for the frequencies shown except in the range from 6 to 12 cps. 

A line with a contracting taper of 6 = -0. 3 differs little from a nontapered line (6 = 0) 
up to 15 cps. For frequencies greater than 15 cps the outlet pressure amplitude pertur- 
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Figure 20. - Variation of resonant frequency with taper factor. 


bation is considerably greater. 

This example indicates that the 
dynamic response of a liquid system 
can be altered significantly with a cer- 
tain amount of control without signifi- 
cantly altering the static response of 
the system; however, the primary 
intent of the example is to illustrate 
the use of the wave -plan analysis. 
Although the example was applied to 
a system with a linear taper and a 
sinusoidal disturbing function, a sys- 
tem of any arbitrary taper with any 
disturbing function could have been 
analyzed with no additional difficulty. 


CONCLUDING REMARKS 

The analytical method presented in this report provides a means of obtaining dis- 
tributed parameter solutions to a variety of unsteady internal flow problems for liquids 
flowing in conduits. An advantage of the method is that the complete solution is obtained. 
For example, both the transient and the steady-state response to a suddenly imposed 
periodic flow disturbance is obtained. Furthermore, the disturbing function can be of 
arbitrary form and need not be periodic. Nonlinear effects are easily included. In addi- 
tion, the wave-plan method is advantageous in making certain types of dynamic response 
calculations. For example, the response of liquid-filled lines having types of axial 
cross-sectional area distributions for which there would be little hope of obtaining 
closed-form analytical solutions can be easily handled. 

The wave-plan analysis can readily solve problems in which there is interaction 
between the structural motion of the conduits and the perturbations in the fluid flowing 
within the conduits. An example of such a problem is the calculation of propellant flow 
in the feed system of a rocket booster experiencing longitudinal oscillations. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, February 4, 1965. 
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APPENDIX A 


COMPUTER SUBROUTINES 

Digital-computer routines have been written to solve for conditions at various line 
discontinuities after wave action in terms of the conditions at the discontinuity prior to 
wave action, the magnitudes of impinging waves, and the physical characteristics of the 
discontinuity. These routines have been formulated in such a manner that they may be 
easily incorporated into a computer program to synthesize different liquid flow systems. 

In all cases the computer routines have been written in nonsubscripted notation, 
which corresponds exactly to the notation presented in the text of this report. In the 
calling vector that calls the subroutine into the main program, the subscripted counter- 
parts of the nonsubscripted subroutine variables are identified. 

This mechanism of identifying subscripted variables in a calling vector makes it 
possible to use the analyses presented in the text (for a contraction or expansion with a 
left-end terminal orifice and flow to the right) to solve for conditions at a diameter dis- 
continuity with a right-end terminal orifice and flow to the left. In this manner the anal- 
yses presented in the text of this report have been generalized. Proper identification of 
the variables in the calling vector is tantamount to orienting the discontinuity to corre- 
spond to the case analyzed in the text. A table incorporated into each subroutine gives 
the proper subscripted-to-nonsubscripted dependence for each case. 

Each of the analyses presented in the text resulted in a polynomial that is quadratic 
in a velocity term. For several reasons these equations were solved by an iterative 
manner by employing Newtonian extrapolation. The roots of the polynomials were deter- 
mined by looping the following equations: 


error = 

2aVll + b 
Vll = Vll - error 

The looping is continued until the error is of sufficiently small magnitude. An ac- 
ceptable solution can be obtained within a very few cycles because the wave-plan solution 
deals with small changes, and the value of the velocity before wave action is a good 
approximation of the final value. In addition, the coefficients b and c are generally 
large and always much larger than a. These conditions lead to a rapidly diminishing 
error term. 

A direct solution of the quadratic equations by the quadratic formula results in a 
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solution that is a small difference between two large numbers. In some cases sufficient 
accuracy cannot be obtained even when double -precision computing techniques are em- 
ployed. A brief description of the program and the Fortran IV program for each com- 
puter routine follows. 


Terminal Orifice Subroutine 


Equations (16) and (17) are the basic equations employed in this analysis. They are 
derived in the text for an orifice on the left end of a line. The table at the beginning of 
the subroutine gives the identities that must be substituted into the calling vector depend- 
ing on whether the orifice is on the right or left end of the line. Before the subroutine 
is entered, the following terms must be identified with a numerical value: 

B orifice coefficient 

HEN reservoir head for left end 
HEX reservoir head for right end 

J time counter 

K number of time increments to nearest discontinuity 

L position counter 

The computer program for a terminal orifice is as follows: 


SUBROUTINE T OR < H2 , DH2 , V2 , Hi 1 , H22 , DH22 , V2 2 ) 
SUBROUTINE TERMINAL ORFICE 
COMMON/BOX/A, B,C, G »BF »ERRDR 

LEFT END RIGHT END 


c 

H2 

HR ( L ♦ J— 1 ) 

HL ( L , J— 1 ) 

c 

V2 

V ( L » J-l ) 

-V ( L , J-l ) 

c 

DH2 

DHL ( L+l » J— K ) 

DHR ( L-l , J 

c 

V22 

V ( L , J ) 

-V ( L , J ) 

c 

H22 

HR (L , J ) 

HL ( L , J ) 

c 

DH22 

DHR ( L , J ) 

DHL ( L , J ) 

c 

HI 1 

HEN 

HEX 


BB = C*B*B/G 

CC=— B*B*( H11-H2-2 .*DH2+C*V2/G ) 


IF(CC) 10,10,4 
4 BB = -BB 
CC = -CC 

10 V2 2 = V 2 

11 ER = (V22*V22+BB*V22+CC)/(2.*V22+BB) 
V2 2 = V22-ER 

IF ( ABS ( ER ) -ERROR ) 12,12,11 


12 DH22 = DH2+C/G* ( V22-V2 ! 
H2 2 =H2 + DH2+DH2 2 
RETURN 
END 
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Friction Orifice Subroutine 


The basic equations for this routine are equations (18) and (19). Before this routine 
is entered, the following terms must have a numerical value: 

BF friction orifice coefficient 

J time counter 

K number of time increments to nearest discontinuity 

L position counter 

The computer routine is the following: 


SUBROUTINE FOR ( Hi »DHl ,DH2 »H2 ,V1 *Hll >DH2 2 »DH11»H22»V11) 
SUBROUTINE FRICTION ORFICE 


c 

HI 

HL ( L , J-l ) 

c 

DH1 

DHR ( L-l ,J-K) 

c 

DH2 

DHL ( L+l , J-K ) 

c 

H2 

HR ( L , J-l ) 

c 

Vl 

V(L , J-l ) 

c 

HI 1 

HL ( L , J ) 

c 

DH11 

DHL ( L , J ) 

c 

DH22 

DHR(L.J) 

c 

H22 

HR ( L , J ) 

c 

Vll 

V ( L , J ) 


COMMON/BOX3/ A, B) C * G * BF »ERROR 
BB = 2.*BF*BE*C/G 

CC=-(H1+2.*DH1+2.*C/G*V1-(H2+2.*DH2) ) *BF*BF 


IF(CC) 3,3.2 

2 BB = -BB 
CC = -CC 

3 VI 1 = Vl 

b ER = ( Vl 1*V1 1+BB*V11+CC)/(3B+2.*V11) 
Vll = Vll-ER 

IF ( ABS < ER ) -ERROR ) 5,5,4 
5 DH11 = DH1+C/G*< Vl-Vll ) 

H 1 1 = H1+DH1+DH1 1 

DH2 2 = DH2+C/G*(V11— Vl ) 

H22 = H2 + DH2 + DH2 2 

RETURN 

END 


Diameter Discontinuity Subroutine 


Equations (20) to (29) are employed in the diameter discontinuity subroutine. These 
equations are derived for flow in the positive direction (to the right) for both an expan- 
sion and a contraction. Before the subroutine is entered, the direction of flow after 
wave action is determined, and this direction in turn determines which subscripted var- 
iables are inserted into the calling vector. 

The subroutine also includes the case where no viscous loss is considered across 
the discontinuity. To utilize this case, it is necessary to assign the term LOSS a value 
of zero. Additional terms that must be defined are the following: 
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J time counter 

K number of time increments to nearest discontinuity 
L position counter 
The computer routine is as follows: 


SUBROUTINE D I SC ( H 1 , Vl » DH1 , H2 , V2 , DH2 , Cl * C2 * A1 » A2 ,H1 1 , DHl 1 *H22 , DH22 , 
1V1 1 »V22 ) 

C DISCONTINUITY SUBROUTINE 

COMMON/BOX/A »B»C* G,BF, ERROR 
C0MM0N/B0X2/ LOSS 

C FLOW TO RIGHT FLOW TO LEFT 

C CCC GREATER THAN 0 CCC LESS THAN 0 


c 


HI 

HL ( A* J-l ) 

HR ( L , J-l ) 

c 


VI 

VL <L,J-1 ) 

-VR ( L » J-l ) 

c 


DHl 

DHR ( L— 1 » J-K ) 

DHL ( L+l , J-K ) 

c 


H2 

HR ( L , J— 1 ) 

HL ( L ♦ J-l ) 

c 


V2 

VR ( L * J— 1 ) 

VL ( L , J— 1 ) 

c 


DH2 

DHL ( L+l * J-K ) 

DHR ( L-l , J-K ) 

c 


HI 1 

HL { L * J ) 

HR ( L , J ) 

c 


Vll 

VL ( L » J ) 

— VR ( L » J ) 

c 


DHl 1 

DHL ( L » J ) 

DHR ( L » U ) 

c 


H22 

HR < L » J ) 

HL(L,J) 

c 


V22 

VR(L.J) 

— VL ( L » J ) 

c 


DH22 

DHR ( L * J ) 

DHL ( L » J ) 

c 


Cl 

CL ( L ) 

CR ( L ) 

c 


C 2 

CR ( L ) 

CL ( L ) 

c 


A 1 

AL ( L ) 

AR (L ) 

c 


A2 

AR ( L ) 

AL ( L ) 



R=A1/A2 




CC= - 

<Hl+2.*DHl+Cl*Vl/G-H2-2 

.*DH2+C2*V2/G) 



BB = 

Cl /G+C 2*R /G 




IF (LOSS) 4 * 3 » 4 



3 

AA = 

(R*R-1. )/(2.*G) 




GO TO 

10 



4 

I F ( R — 

1 . ) 1,1.2 



1 

AA = 

( R*R-R ) /G 




GO TO 

10 



2 

AA = 

( 3.*R*R-2.*R~1 , ) / (4.*G) 



10 

V 1 1 = 

Vl 



5 

ER = 

( AA*V11*V11+BB*V11+CC ) /( 2.*AA*V11+BB ) 



VI 1 =V 1 1-ER 



IF (ABS(ER)— ERROR) 6»6>5 
6 DHl 1 = DH1+C1 /G*(V1-V11 ) 
V22=R*V11 

DH22 = DH2+C2 /G*(V22-V2 ) 

HI 1 =H1+DH1+DH 1 1 

H22=H2+DH2+DH22 

RETURN 

END 


33 



APPENDIX B 


DIGITAL MODEL OF HYDRAULIC LINE WITH OSCILLATING INPUT ORIFICE 

A schematic drawing of the system chosen for this example is shown in figure 11 

(p. 20). 

The system consists of a variable-area orifice at A that perturbs the system and a 
fixed orifice at B. The conduit is bounded by pressure reservoirs. (The reservoir 
pressures can vary. ) 

To keep the solution general, N - 2 friction orifices are inserted, which result in 
N discontinuities. A time interval is so chosen that each discontinuity is one time in- 
terval from adjacent discontinuities. Thus 

At* = k (-L- 
C \N - 1 

This time increment is a multiple of the working time increment At, which is 
chosen small enough so that all disturbing functions can be accurately described by step 
functions. 

Initial values for velocity through each discontinuity and the pressure to the right 
and left of each discontinuity must be inserted into the computer as data or the computa- 
tion of these quantities from the initial steady condition must be integrated into the com- 
puter program. 

The disturbing function can either be read into the computer as a function of time or 
it can be computed. For this case a sinusoidal orifice coefficient for the inlet (at A) is 
taken to be 


B = BO + BA SIN(27r FJAt) 

The orifice coefficient for the outlet is considered to be constant and is denoted 
by BN. 

The following data are needed as input to the computer: 

C wave velocity in conduit, ft/sec 

L length of line, ft 

BO mean orifice coefficient when subjected to sinusoidal perturbations 
BOl steady-state inlet orifice coefficient 

BA amplitude of inlet orifice coefficient perturbation 

F frequency of pressure perturbation 
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BN outlet orifice coefficient 

VO steady line velocity, ft/sec 

HEN pressure head of inlet reservoir, ft 
HEX pressure head of outlet reservoir, ft 

g acceleration due to gravity, ft/sec 

N number of discontinuities 

NN number of working time increments for wave to travel length of conduit 

K number of working time increments between discontinuities 

M total number of time increments for which computations will be made 

N2 number of friction orifices for which pressure and velocity data will be 

printed 

In setting up the computer program the notation j = x,y, z is used. The notation 
means that computations are carried out for j starting at j = x. The computations are 
repeated for values of j that are increased by y until j > z. 

The following output parameters are printed out by the computer: 

V(1,J), HR(1, J) 

V(N, J), HL(N, J) 

V(N2, J), HR(N2, J), HL(N2, J) 

For a specific example, the fluid system constants appearing in the computer pro- 
gram were chosen nominally to match experimental conditions in the Lewis line- 
dynamics facility shown schematically in figure 11(b) (p. 20). Two sets of line condi- 
tions were used. The first (case I) 

table n. - constants for example was chosen in conjunction with an 

Case I 

3800 
68 
0.8 

0.48,0.32,0. 20,0. 15 
1.0 

14 
520 
0 

32. 2 
17 

1 

170 

16 

40 

5 

0. 00001 


C, ft/sec 
L, ft 

BOl = BO 

BA 

BN 

VO, ft/sec 
HEN, ft 
HEX 

g, ft/sec 1 2 * * 5 
N 

K 

M 

NN 

F, cps 
N2 

ERROR 


Case n 


0.65 

0. 08,0.22,0.36,0.50 
0. 42 


analytical study of the effect of the 
size of the input amplitude on the 
linearity of the perturbation imposed 
on the mean flow and pressure in the 
line. The second (case II) was chosen 
to match conditions frequently obtained 
in the line. 

The constants listed in table H 
were employed. The Fortran IV com- 
puter program is the following: 
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COMMON/BOX/A, S.C. G.SF, ERROR 

DIMENSION V( 20,101 ) ♦DHL t 20. HO) * DHR ( 20, 110). HR (20,101). HL( 20*101 ) » 
1 B 2 < 101 ) *X (550) > Y1 (550) » Y2 ( 550 ) > Y 3 ( 550 ) * Y4 ( 550 ) » Y5 ( 550 ) » Y6 ( 550 ) >Y7 ( 
1550) * Y 8 t 5 5 0 ) . Y 9 ( 5 5 G ) » Y 1 0 ( 550) > Y 1 1 (550) 

G = 32 . 2 

KKK = 1 0 1 

C 
C 
C 
C 
C 
C 
C 
C 
C 

c 
c 
c 


I- UK TRAN 

HL(L.J-l) J = 2 » M+l 
HR { L » J — 1 ) 

DHR ( L- 1 * J- 1 ) 

DHL ( L+l * J— 1 ) 

V ( L » J-l ) 

HL ( L * J ) 

HL(L.J) 

HR ( L , J ) 

V ( L * J ) 

DHR ( L . JK ) 

DHL (L.JK) 


ri l wni_ 




HL(L.J-l) J=1.M 
HR ( L » J-l ) 

DHR ( L-l * J-K ) 
DHL! L+l , J-K) 

VI L » J-l ) 

V ( L * J ) 

HL ( L » J ) 

HR ( L . J ) 

V ( L * J ) 

DHR < L , J > 

DHL t L » J ) 


C READ PROGRAM CONSTANTS 

1 READ! 5 » 100)C»XL.B0»SA.F>3N,V0 ♦HEN .HEX .BDl .ERROR 
READ! 5» 101 )N»NN»K »M,N2 

100 FORMAT(8F10.5) 

101 FORMAT! 1015) 

C WRITE PROGRAM CONSTANTS 

WRITE (6. 901) 

WRITE(6*900)C»XL .80 , 3A , F , BN , VO .HEN , HEX . BD 1 

900 FORMAT ( 12G10. 3 ) 

901 FORMAT ( 1H1 .5X.2HC-,5X,2HL=»8X,3HBD=,7X»3HBA=»7X»2HF=,3X»3HBN=,7X, 
1 3HV0=,7X.AHHEN=.6X»AHHEX=.6X ,AHB01= ) 

WRITEI6.902) 

9C2 FORMAT! 3X , 2HN= , 2 X > 3HNN= . 3 X , 2HK= , 3X , 2HM = , 2X , 3HN2= ) 

WR I T E ( 6 » 1 0 1 ) N . NN . K , M . N2 

C COMPUTE INITIAL CONDITIONS 

L2 = -l 
XN = N 
DO 5 L=1»N 

V t L » 1 ) = VO 
DO 5 J= 1 , K 
DHL ( L » J ) = • 0 

5 DHR ( L , J ) = .0 

HR ( 1 ,1 ) =HEN- ( V0/B01 ) **2 
HL(N.1)=HEX+(V0/BN)**2 

B F = VO / ( (HR(l.l)-HL(N.l) )/( XN-2 .))**. 5 
CONST = VO*VC/ ( BF + BF ) 

N 3 = N - 1 

DO 6 L=2,N3 

HRfL.l )=HR(L-1 ,1 ) -CONST 

6 HL ( L .1 ) =HR ( L-l .1 ) 

X N N = N N 

DT=XL/ ( XNN*C ) 

NK=M +1 
N K 2 = 2*NK 
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WRITE INITIAL CONDITIONS 


WRITE (6,102) 

102 FORMAT (1H0»5X»8HHL(L*0)=) 

WRITE (6,103) (L,HL(L,1 ) »L=2,N> 

103 FORMAT !4( I6*F16.4,5X) ) 

WRITE (6,104) 

104 FORMAT ( 1HO»5X,8HHR(L»0)= ) 

WRITE (6,103) (L,HR(L,1) ,L = 1,N3! 

WRITE (6,105) BF 

105 FORMAT ( 1H0»5X»3HBF=,E14,6 ) 

MN = M 

MN1=MN+1 

C0N=2.0*3.141593*F*DT 

Xl=V(l,l) 

X2 = HR ( 1,1) 

X4 = HL( N2, 1 ) 

X5 = HR ( N2 , 1 ) 

X6= V ( N ,1) 

X7 = HL ( N , 1 ) 

X8=B0+BA*SIN (CON ) 

C START TIME LOOP 

20 DO 7 J = 2 , MN 1 
AJ = J+L2 

C COMPUTE CONDITIONS AT LEFT END 

B2(J) = BO+BA*SIN(CON*AJ) 

B=B2 ( J ) 

JK = J + K-l 

CALL TOR ( HR ( 1 , J-l ) , DHL ( 2 , J- 1 ) , V ( 1 , J- 1 ) , HEN, HR ( 1 , J ) ,DHR ( 1 , JK ) , 

1 V ( 1 , J ) ) 

C COMPUTE CONDITIONS AT FRICTION DRIFICES 

DO 8 L = 2 , N 3 

8 CALL F0R(HL(L,J-1 ) ,DHR(L-1,J-1 ) ,DHL!L + 1 *J-1 ) ,HR(L,J-1 ) ,V( L , J-l ) , 
1HL ( L » J ) »DHR ( L , JK ) »DHL ( L , JK) »HR ( L , J ) *V ( L , J) ) 

V2 2 =-V ( N, J-l ) 

C COMPUTE CONDITIONS AT RIGHT END 

B = BN 

CALL TOR ( HL ( N , J-l ) ,DHR ( N-l , J-l ) , V22 , HEX , HL ( N , J ) , DHL ( N , JK ) , 

1 V ( N , J ) ) 

V ( N , J ) =-V ( N , J ) 

IF ( J-KKK ) 7,9,9 
7 CONTINUE 
GO TO 10 

C STORE RESULTS FOR PRINTING 

9 DO 50 J = 2 , KKK 
LL= J+L2 

Y 1 ( LL ) = V(1,J) 

Y2(LL)= HR ( 1 , J ) 

Y4( LL) = HL ( N 2 , J ) 

Y5 ( LL ) = HR ( N 2 , J ) 

Y6 ( LL) = V ( N , J ) 
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Y8 ( LL ) = B2 ( J ) 

50 Y7 ( LL ) = HL ( N » J ) 

L2=L2+KKK -1 

200 FORMAT ( 7E16. 8 ) 

C INITIALIZE 

DO 15 L=1,N 

V ( L » 1 ) = V ( L * KKK ) 

DO 15 J= 1 » K 
KK =KKK+J-1 
DHL(L»J)=DHL(L»KK) 

15 DHR(L,J)=DHR(L,KK) 

HR ( 1 .1 ) =HR ( 1 » KICK ) 

HL(N»1 ) =HL ( N * KKK ) 

DO 16 L = 2 * N 3 

HR ( L » 1 ) =HR ( L » KKK ) 

16 HL ( L > 1 ) =HL ( L * KKK ) 

MN = MN-KKK+1 

MN 1 =MN+1 
IF(MN) 19,19,20 

C STORE RESULTS FOR PRINTING 

10 DO 55 J = 2 * MN 1 
LL= J+L2 

Y1 ( LL ) = V ( 1 , J ) 

Y2(LL)= HR ( 1 ,J) 

Y A ( L L ) = HL ( N2 , J ) 

Y 5 ( LL ) = HR ( N2 , J ) 

Y7 ( LL ) =HL ( N , J ) 

Y6(LL)=V(N,J) 

55 Y 8 ( LL) =B2 ( J) 

19 NK=NK-1 

WRITE (6,201) 

201 FORMAT ( 1 HI , 5 X , 2 HV 1 , 1 6 X , 3 HV-N , 1 3 X , 5 HHL-N2 , 1 IX, 4HHL-N, ,12X,1HB,15X, 
14HHR-1 , 8X » 5HHR-N2 ) 

V.'RI TZ( 6, 111 ) ( Y1 ( J ) ,Y6 ( J ) ,YA ( J ) ,Y7 ( J) ,Y8 ( J ) ,Y2 ( J) ,Y5( J) , J = 1 ,NK) 

1 1 1 FORMAT ( 7F16.4 ) 

DO 56 1 = 1, NK 

Y 1 ( I ) =Y1 ( I ) — X 1 
Y2 ( I ) = Y2 ( I ) -X 2 
Y 4 ( I ) = Y4 ( I ) — X4 
Y5 ( I ) = Y 5 ( I ) -X 5 
Y 6 ( I ) =Y6 ( I ) — X6 
Y8 ( I ) =Y8 ( I ) -X8 

56 Y7 ( I ) =Y7 ( I ) -X 7 

C WRITF OUTPUT 


WRITE (6,201) 

WRITE(6,111)(Y1(J),Y6(J),Y4(J),Y7(J),YR(J),Y2(J),Y5(J),J=1,NK) 

GO TO 1 

END 
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APPENDIX C 


DIGITAL MODEL FOR TAPERED HYDRAULIC LINE 

The system chosen for the tapered hydraulic line is shown in figure 17 (p. 26). This 
system consists of a variable -area orifice at A and a fixed orifice at B connected by a 
tapered line with a linear diameter variation between A and B. The line is bounded by 
pressure reservoirs. The tapered line was approximated by a line composed of a dis- 
crete number of constant-diameter sections of decreasing (or increasing) diameter. 

To keep the solution general, N - 2 diameter discontinuities are inserted, which 
result in N discontinuities and N - 1 sections of straight conduit. 

The conduit is so divided that the wave travel times between adjacent discontinuities 
are multiples of the working time interval, and the diameter of the approximating 
straight section of conduit is equal to the average diameter of the section that it is 
approximating. 

It is assumed that the wave velocity remains constant over the entire conduit, which 
is equivalent to assuming that the ratio of the conduit diameter to wall thickness remains 
constant (a reasonable assumption from a structural viewpoint). The conduit is thus 
approximated by N - 1 sections of equal length. The diameter of the approximating 
straight section is given by 


/DB - DaV2I - l\ 

D<R = DA + (— r — )(— ) 

Initial values for velocity and pressure to the right and left of each discontinuity 
must be inserted into the computer, or the equations for the computation of these quan- 
tities from initial steady conditions must be integrated into the computer program. 

The disturbing function can either be read into the computer as a function of time 
or it can be computed. For this model a sinusoidal orifice coefficient for the inlet (at A) 
is taken to be 


B = BO + BA SIN(2irFJAt) 

The orifice coefficient for the outlet is constant and denoted by BN. 
The following data are needed as input to the computer: 

DA diameter at A, ft 
DB diameter at B, ft 
CA wave velocity in conduit, ft/sec 
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L length of line, ft 

BOl steady-state input orifice coefficient 

BO mean orifice coefficient when subjected to sinusoidal perturbations 
BA amplitude of orifice coefficient perturbation 

F frequency pf pressure perturbation 

BN output orifice coefficient 

HEN pressure of input reservoir, ft 

QO steady-state flow rate, ft 

g acceleration due to gravity, ft/sec 

N number of discontinuities 

NN number of working time increments for wave to travel length of conduit 
K number of working time increments between discontinuities 

M total number of time increments for which computations will be made 
N2 number of diameter discontinuity for which pressure and velocity data will 
be printed 

The following output parameters are printed out by the computer: 

VR( 1, J) , HR(1, J) 

VL(N, J), HL(N, J) 

VL(N2,J), HL(N2,J) 

HR(N2, J) 

For a specific example the effect of line taper on the response of a particular sys- 
tem of the type shown in figure 17 (p. 26) was computed. Calculations were made for 
nine different degrees of line taper. The specific data used in this analysis were as 


follows: 


TABLE m. - END DIAMETERS AND 


CONDITIONS FOR TAPER STUDIES 


DA 

DB 

BO 

BA 

BN 

1 

1 

1 

0. 1 

0.2 

.9 

1. 1 

1.235 

. 1235 

. 165 

.8 

1. 2 

1. 562 

. 1562 

. 139 

. 7 

1.3 

2. 041 

. 2041 

. 1185 

.6 

1.4 

2.78 

.278 

. 102 

. 5 

1.5 

4. 0 

. 4 

.089 

1. 1 

.9 

.825 

.0825 

. 247 

1.2 

.8 

.695 

.0695 

.313 

1.3 

.7 

. 593 

.0593 

.408 


CA = 3000 ft/sec 
L = 50 ft 
HEN = 660 ft 
QO = 6. 4 cu ft/sec 
N = 11 
NN = 10 
M = 300 

Table HI gives the end diameters and conditions 
for the various tapers studied. 

These cases were studied for frequencies of 
5, 7.5, 10, 12.5, 15, 17.5, 20, 22.5, and 25 cps. 

The Fortran IV program for this example is as 
follows: 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

COMMON/BOX/A,B,C» G,BF, ERROR 
C0MM0N/B0X2/ LOSS 

DIMENSION D { 1 00 ) »VR(20,101) »HR(20*101) *CL(XOO) »CR(100)»VL(20»201) » 
1 Y2 ( 500 ) , Y3 ( 5 00 ) »Y4(500),Y5(500),Y6(500) » 

2Y1(500),AL(500),AR(500),HL(20,101!»DHL(20»110) ,DHR(20,110) 

PI = 3.141593 
G = 3 2 • 2 


FORTRAN 
HR ( L » J— 1 ) 
DHL ( L ,J-1 ) 

V ( L » J — 1 ) 
HR (L.J) 
DHR(L,JK) 
V( L.J ) 
HL(L.J-l) 
HL (L.J) 
DHL (L.JK) 


AL I UnL I'lCHIVll'IU 


HR ( L , J-l ) 
DHL ( L . J-K) 
VtL.J-l ) 

HR ( L . J ) 
DHR(L.J) 

V ( L » J ) 

HL ( L , J-l ) 
HL ( L . J ) 
DHL (L.J) 


J = 1 »M 


C READ DATA 

600 READ ( 5. 100 ) DA.DB.BOl »B0 , B A , BN . HEN , XL .C » F , L OSS » ERR OR , QO 
100 FORMAT ( 8F10. 5 ) 

READ ( 5, 101 1 K.M.N ,N2 ,NN 


C WRITE INPUT DATA 

WR ITF(6»113) 

113 FORMAT ( 1H1 , 2 X . 3HDA= , 8 X . 3 HDB = , 9X . 4HB0 1 = , 8X , 3HB 0= , 9X ,-3HB A = ,9X»3HBN=, 
19X.4HHEN=»8X.3HXL=.9X.3HC= .9X,2HF=»9X,3HQ0=) 

WR I TE ( 6, 112) DA.DB.BOl, BO, BA, BN, HEN, XL, C ,F,00 
112 FORMAT ( 11F11 .4) 

101 FORMAT (815) 

WR I TE ( 6 , 1 16 ) 

116 FORMAT ( 1H0,5X,2HK=,6X,2HM=,8X,2HN=,8X,3HN2=»7X»3HNN= ) 

WR I T E ( 6 , 1 1 5 ) K,M,N,N2,NN 
115 FORMAT (5110) 

MN = M 
L 2 = 0 

N 1 = N-l 
XNN=NN 

DT = XL /C /XNN 
C0N=2.*PI*F 
XN = N 
KKK = 101 


C 


COMPUTE INITIAL CONDITIONS AND LINE CONSTANTS 


DO 5 1=1 , N 1 
A I = I 

5 D ( I ) = DA+(DB-DA)*(2.*AI-1. )/(2.*(XN-l. ) ) 
VL ( 1 , 1 ) = QD*4./(PI*D(1)*D(1> ) 

HR (1,1) = HEN- (VL ( 1 ,1 ) /B31 ) **2 
DO 6 L=2 ,N1 
CL ( L ) = C 
CR ( L ) = C 

AL ( L ) = PI*D(L-1 )*-»2*.25 
AR ( L ) = PI*D(L )**2*.25 
VL ( L, 1 ) = OO/AL ( L ) 

VR ( L , 1 ) = QO/AR(L ) 

R = AL ( L ) /AR ( L ) 
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IF (R-l.) 8,8,3 
8 AA = (R*R-R)/G 
GO TO 4 

3 AA = (3.*R*R-2.*R-1.)/<4.*G) 

4 HL ( L , 1 ) = HR ( L — 1 , 1 ) 

6 HR l L , 1 ) = HL(L,1 )-AA*VL(L,l )**2 
VR ( N , 1 ) = QO/AR ( N-l ) 

HL ( N , 1 ) = HR ( N-l , 1 ) 

HEX = HL( N,1 )-( VR (N ,1 ) /8N)**2 
MN = M+l 
DO 7 L = 1 » N 
DO 7 J= 1 , K 
DHL ( L , J ) = .0 

7 DHR (L, J) = .0 
WRITE (6,107) HEX 

107 FORMAT ( 2H0V , 1 6X , 2HHL , 1 6X , 2HHR , 5X , 4HHEX = , F 1 2 . 4 ) 

108 FORMAT ( 3E16. 8 , I 5 ) 

X 1 = VL ( 1,1) 

X2 = VR ( N ,1) 

C WRITE INITIAL CONDITIONS 

WRITE (6,108) (VL ( L— 1 ,1 ) , HL ( L , 1 ) »HR(L-1»1) » L » L = 2 »N ) 

X3 = HR (1,1) 

X4=HL ( N ,1) 

X5 = HL ( N2 , 1 ) 

X6 = HR ( N2 , 1 ) 

C SET UP TIME LOOP 

30 DO 9 J=2 ,MN 
JK = J + K-l 
A J= J+L2-1 

C DEFINE CONDITIONS FOR LEFT END 

B = BO + BA*SIN(CON*DT*AJ) 

L = 1 

C TERMINAL ORFICE SUBROUTINE 

CALL TOR ( HR ( L , J-l ) ,DHL ( L + l, J-l ) ,VL ( L, J-l ) , HEN ,HR ( L , J ) , DHR ( L, JK) , 

1 VL ( L, J ) ) 

DO 10 L = 2 » N 1 

C CHECK CONDITIONS FOR DIAMETER DISCONTINUITY 

CCC =HR ( L , J-l ) +2 . * DHL ( L-l , J-l ) + CL ( L ) *VL (L , J-l ) /G-HL (L , J-l )-2.* 
1DHR (L + l , J-l ) +CR ( L ) *VR < L » J-l ) /G 
IF(CCC) 1,1,2 
1 V2=-VL(L, J-l ) 

V1=-VR(L, J-l ) 

C DIAMETER DISCONTINUITY SUBROUTINF 

CALL DISC(HR(L,J-1) ,V1 , DHL ( L+l , J-l ) , HL ( L , J-l ) , V2 , 

1DHR ( L-l , J-l ) ,CR(L),CL(L),AR(L),AL(L),HR(L,J) ,DHR(L,JK) ,HL!L»J) , 
2DHL ( L , JK ) , VR(L,J>* VL(L,J!) 

VR ( L , J ) =- VR ( L , J ) 

VL ( L , J ) = — VL ( L , J ) 

GO TO 10 
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2 CALL DISC ( HL ( L , J-l ) ,VL ( L, J-l ) ,DHR ( L-l , J-l ) .HR (L , J-l ) ,VR( L, J-l ) > DHL 

1 (L+l » J-l ) ,CL ( L ) ,CR(L > ,AL ( L) ,AR(L ) ,HL(L ,J) .DHL f L , JK) ,HR(L , J ) »DHR ( L , 
2JK) ,VL(L,J) ,VR(L»J) ) 

10 CONTINUE 

C DEFINE CONDITIONS FOR RIGHT END 

L = N 
B=BN 

V2 = -VR ( L * J-l ) 

C TERMINAL ORFICE SUBROUTINE 

CALL TOR <HL< L .J-l > »DHR (L-l» J-l ) ♦ V2 . H EX . HL ( L , J ) . DHL ( L . JK ) , 

1 VR t L . J ) ) 

VR(L»J)=— VR(L.J) 

IF(J-KKK) 9.26.26 
9 CONTINUE 

C END TIME LOOP 


C STORE RESULT FOR PRINTING 

26 DO 27 J = 1 » KKK 
LL = J+L 2 

Y 1 ( LL ) = VL ( 1 , J ) 

Y2 ( LL ) =VR ( N .J ) 

Y3 ( LL ) =HR ( 1 , J ) 

YA ( LL ) =HL ( N .J ) 

Y5 ( LL ) =HL ( N2 , J ) 

27 Y6 ( LL ) =HR ( N2 » J ) 

L2=L2+KKK-1 

C INITIALIZE 

DO AO L = 1 » N 

VL ( L , 1 ) = VL ( L .KKK ) 

VRCL.l ) =VR ( L , KKK ) 

DO AO J = 1 , K 

KK = KKK + J- 1 

DHL ( L , J ) = DHL ( L , KK ) 

AO DHR(L,J)=DHR(L,KK) 

DO 28 L = 1 . N 
HR(L.1)=HR(L.KKK) 

28 HL(L.1)=HL(L.KKK) 

MN=MN-KKK+ 1 
IF(MN) 29.29,30 

29 NK=M+ 1 

C WRITE OUTPUT 

WRITE (6,109) 

109 FORMAT ( 3H0V1 »16X .3HHR1 , 1AX.AHVN » 1 2X » 5HHLN .11X.AHHLN2.12X, A HHRN 
12 ) 

WRITE(6»110) ( Y 1 ( J ) » Y 3 ( J ) ,Y2(J),YA(J) » Y 5 ( J ) ,Y6(J) , J = 1 » NK ) 

DO 31 J = 1 » NK 
Y 1 ( J ) = Y1 ( J ) — X 1 
Y2( J)=Y2( J1-X2 
Y3 ( J ) =Y3 ( J ) — X3 
YA ( J ) =YA ( J ) -XA 
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Y5( J)=Y5( J)-X5 
31 Y6( J)=Y6( J)-X6 
WRITE (6,109) 

WRITE(6,110)(Y1(J),Y3(J),Y2(J),Y4(J),Y5(J),Y6(J),J=1,NK) 
110 FORMAT ( 6G16. 8 ) 

GO TO 600 
END 
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APPENDIX D 


SYMBOLS 


A 

line area, sq ft 

E c 

elastic modulus of conveyer, 

A(i) 

line area of i, sq ft 


lb/sq ft 

A(j) 

line area of j, sq ft 

E f 

elastic modulus of fluid, lb/sq ft 

A o 

orifice area, sq ft 

f 

Darcy friction factor 

A 1 

line area to left of diameter dis- 

/ 

function 


continuity, sq ft 

g 

acceleration due to gravity, 

A 2 

line area to right of diameter 


ft/sec 2 


discontinuity, sq ft 

H 

pressure head, ft 

a, b, c 

polynomial coefficients 

HEN 

pressure head of inlet reservoir, 

B 

orifice coefficient, (ft/sec^)^/ 2 


ft 

b f 

friction orifice coefficient, 
(ft/sec^)'*’/ 2 

HEX 

pressure head of outlet reser- 
voir, ft 

B 1 

orifice coefficient before wave 
action, (ft/sec 2 ) 1//2 

H 1 

pressure head to left of discon- 
tinuity before wave action, ft 

B 2 

orifice coefficient after wave 
action, (ft/sec 2) ^ 2 

H 2 

pressure head to right of discon- 
tinuity before wave action, ft 

C 

C D 

sonic velocity in line, ft/sec 
orifice discharge coefficient 

H 11 

H 22 

pressure head to left of discon- 
tinuity after wave action, ft 

pressure head to right of discon- 

C(i) 

wave velocity in line i, ft/sec 

tinuity after wave action, ft 

C(j) 

wave velocity in line j, ft/sec 

AH 

pressure head change across 

C 1 

wave velocity to left of discon- 


wave or orifice, ft 


tinuity, ft/sec 

AHj 

wave impinging from left of dis- 

C 2 

wave velocity to right of discon- 
tinuity, ft/sec 


continuity before wave action, 
ft 

D 

line diameter, ft 

ah 2 

wave impinging from right of 

D 

mean diameter of tapered line, 
ft 


discontinuity before wave 
action, ft 
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AH 


11 


AH 


22 


Ah L 

J 

K 

L 

P 

AP 

Q 

R 

R(i) 

T(i) 

t 

At 

t" 

t + 


wave leaving discontinuity on 
left side after wave action, ft 

wave leaving discontinuity on 
right side after wave action, 
ft 

pressure head loss over small 
line length, ft 

time subscript 

number of working time incre- 
ments between discontinuities 

line length, ft 

pressure in line, Ib/sq ft 

pressure wave, lb/sq ft 

volume flow rate, cu ft/sec 

area ratio at diameter discon- 
tinuity 

reflection coefficient in line i 
transmission coefficient in line i 
time, sec 

incremental time interval, sec 

short time before wave action, 
sec 

short time after wave action, sec 


V 

velocity in line, ft/sec 

VE X 

average velocity of moving end 
before wave action, ft/sec 

ve 2 

average velocity of moving end 
after wave action, ft/sec 

VO 

velocity in line adjacent to mov- 
ing orifice, relative to orifice, 
ft/ sec 

V 1 

velocity in line to left of discon- 
tinuity before wave action, 
ft/sec 

V 2 

velocity in line to right of discon- 
tinuity before wave action, 
ft/sec 

V 11 

velocity in line to left of discon- 
tinuity after wave action, 
ft/ sec 

V 22 

velocity in line to right of discon- 
tinuity after wave action, 
ft/sec 

X 

position, ft 

6 

line taper factor (eqs. (31)), ft 

P 

g 

mass density of fluid, slugs/ft 

T 

wall thickness of conduit, ft 
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